%%%-------------------------------------------------------------------
%%% File    : p23.erl
%%% Author  : Plamen Dragozov <plamen at dragozov.com>
%%% Description : 
%%% A perfect number is a number for which the sum of its proper 
%%% divisors is exactly equal to the number. For example, the sum of the 
%%% proper divisors of 28 would be 1 + 2 + 4 + 7 + 14 = 28, which means 
%%% that 28 is a perfect number.
%%% 
%%% A number whose proper divisors are less than the number is called 
%%% deficient and a number whose proper divisors exceed the number is 
%%% called abundant.
%%% 
%%% As 12 is the smallest abundant number, 1 + 2 + 3 + 4 + 6 = 16, the 
%%% smallest number that can be written as the sum of two abundant numbers 
%%% is 24. By mathematical analysis, it can be shown that all integers 
%%% greater than 28123 can be written as the sum of two abundant numbers. 
%%% However, this upper limit cannot be reduced any further by analysis 
%%% even though it is known that the greatest number that cannot be 
%%% expressed as the sum of two abundant numbers is less than this limit.
%%% 
%%% Find the sum of all the positive integers which cannot be written as 
%%% the sum of two abundant numbers.
%%% 
%%% Created : 10 Dec 2008
%%%-------------------------------------------------------------------
-module(p23).

%% API
-compile(export_all).

%%====================================================================
%% API
%%====================================================================

% A list of abundand numbers:
% Abundand = (A1, A2, A3, ..., AN)
% All numbers that are a sum of 2 abundand:
% SumsOfAbundand = set(2A1, 2A2, 2A3, ...., 2AN, (A1+A2),..(A1+AN), (A2 + A3), ..., (A2+AN), ...(AN-1 + AN))
% The sum of all numbers that are not a sum of 2 abundand
% SumOfSums = Sum(1, Max) - Sum(SumsOfAbundand)
% i.e. The sum of all numbers up to Max minus the sum of all numbers
% that are a sum of 2 abundant.

%%--------------------------------------------------------------------
%% Function: solution() ->
%% Description: Returns the sum of all numbers that are NOT a sum of 2
%% abundant numbers.
%%--------------------------------------------------------------------
solution() ->
    sum(28123).

%%====================================================================
%% Internal functions
%%====================================================================
sum(N) ->
    Sum = (N*(N+1)) div 2,
    scan(N, N, ets:new(lookup, []), [], Sum ).

scan(0, _, _, _, Sum)->
    Sum;
scan(N, Max, LookupTable, Acc, Sum) ->
    case is_abundand(N) of
        true ->
            NewAcc = [N|Acc],
            NewSum = update_sum(Sum, N, Max, LookupTable, NewAcc), 
                                                %find all new numbers that are equal to the sum of N and
                                                %another abundand number and subtract them from the sum. 
            scan(N - 1, Max, LookupTable, NewAcc, NewSum);
        _ -> scan(N - 1, Max, LookupTable, Acc, Sum)
    end.

update_sum(Sum, N, Max, LookupTable, [H|T]) ->
    Val = N + H,
    case Val > Max of 
        true -> Sum;%stop the search
        _ ->
            Old = ets:lookup(LookupTable, Val),
            case Old of
                [] ->
                    ets:insert(LookupTable, {Val, true}),
                    update_sum(Sum - Val, N, Max, LookupTable, T);
                _ -> update_sum(Sum, N, Max, LookupTable, T)
            end
    end.

%Check if a number is abundand            
is_abundand(N) ->
    Sum = divisors_sum(N),
    Sum > N.

divisors_sum(N) ->
    {D, P} = divisor_powers(N),
    divisors_sum(D, P, 1) - N.

pow_sum(X, N) ->
    pow_sum(X, 1, N, 1, 1).

divisors_sum([], [], Acc) ->
    Acc;
divisors_sum([DivsH| DivsT], [PowersH| PowersT], Acc) ->
    divisors_sum(DivsT, PowersT, Acc*pow_sum(DivsH, PowersH)).

pow_sum(X, Max, Max, Exp, Acc) ->
    Acc+ X*Exp;
pow_sum(X, Pow, Max, Exp, Acc) ->
    Next = Exp * X,
    pow_sum(X, Pow + 1, Max, Next, Acc + Next).

divisor_powers(N) ->
    {D, P} = divisor_powers(N, 2, 0, [], []),
    {lists:reverse(D), lists:reverse(P)}.

divisor_powers(N, I, Counter, Acc1, Acc2) when I > (N div 2)->
    case true of
        true when I =:= N -> {[N|Acc1], [Counter + 1 | Acc2]};
        true when Counter > 0 -> {[N|[I|Acc1]], [1|[Counter|Acc2]]};
        _ -> {[N|Acc1], [1|Acc2]}
    end;
divisor_powers(N, I, Counter, Acc1, Acc2) ->
    case true of
        true when (N rem I) =:= 0 ->
            divisor_powers(N div I, I, Counter + 1, Acc1, Acc2);
        true when Counter > 0 ->
            divisor_powers(N, I + 1, 0, [I|Acc1], [Counter|Acc2]);
        _ ->
            divisor_powers(N, I + 1, 0, Acc1, Acc2)
    end.
